<!DOCTYPE html>
<html lang="en">

<head>
    <title>QSKJ_LOCAL</title>
    <script src="/Public/js/tianditu.api.js"></script>
    <script src="/Public/js/jquery-3.7.1.min.js"></script>
    <script src="/Public/js/geoTools.js"></script>
    <link rel="stylesheet" href="/Public/css/1.css">
</head>

<body>

    <script>

        function gaussProjectionArcMap(B, L) {
            // CGCS2000 椭球参数
            const a = 6378137.0;          // 长半轴
            const f = 1 / 298.257222101;  // 扁率
            const b = a * (1 - f);        // 短半轴
            const e2 = (a * a - b * b) / (a * a); // 第一偏心率的平方
            const e12 = (a * a - b * b) / (b * b); // 第二偏心率的平方

            // 计算 3度带中央子午线 L0
            const zone = Math.floor(L / 3 + 0.5);
            const L0 = zone * 3;

            // 角度转弧度
            const rad = Math.PI / 180;
            const B_rad = B * rad;
            const L_rad = L * rad;
            const L0_rad = L0 * rad;

            // 经差 l = L - L0（弧度）
            const l = L_rad - L0_rad;

            // 计算辅助量
            const sinB = Math.sin(B_rad);
            const cosB = Math.cos(B_rad);
            const t = Math.tan(B_rad);
            const W = Math.sqrt(1 - e2 * sinB * sinB);
            const N = a / W; // 卯酉圈曲率半径
            const eta2 = e12 * cosB * cosB;

            // 计算子午线弧长 X（幂级数展开）
            const A0 = 1 + (3 * e2 / 4) + (45 * e2 * e2 / 64) + (175 * e2 * e2 * e2 / 256);
            const A2 = (3 * e2 / 8) + (15 * e2 * e2 / 32) + (525 * e2 * e2 * e2 / 1024);
            const A4 = (15 * e2 * e2 / 256) + (105 * e2 * e2 * e2 / 1024);
            const A6 = 35 * e2 * e2 * e2 / 3072;

            const X = a * (1 - e2) * (
                A0 * B_rad - A2 * Math.sin(2 * B_rad) + A4 * Math.sin(4 * B_rad) - A6 * Math.sin(6 * B_rad)
            );

            // 计算 x, y
            const l2 = l * l;
            const l4 = l2 * l2;
            const l6 = l4 * l2;

            const x = X +
                (N * t / 2) * l2 +
                (N * t / 24) * (5 - t * t + 9 * eta2 + 4 * eta2 * eta2) * l4 +
                (N * t / 720) * (61 - 58 * t * t + t * t * t * t) * l6;

            const y =
                N * l +
                (N / 6) * (1 - t * t + eta2) * l * l2 +
                (N / 120) * (5 - 18 * t * t + t * t * t * t + 14 * eta2 - 58 * eta2 * t * t) * l * l4;

            // Y坐标加常数 500,000 并加带号（如19带 → 19,500,000）
            const y_final = y + 500000 + zone * 1000000;

            return {
                x: parseFloat(x.toFixed(6)),
                y: parseFloat(y_final.toFixed(6)),
                zone: zone
            };
        }

        // 测试：北京天安门 (39.9075°N, 116.3972°E)
        const result = gaussProjectionArcMap(31.145,102.549);
        console.log(result); // 应输出 {x: 3447180.30359, y: 19552352.711796, zone: 19}
    </script>
</body>

</html>